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In the theory of Bethe-ansatz integrable qnantum systems, rapidities play an important role 
as they are used to specify many-body states, apart from phases. The physical interpretation of 
rapidities going back to Sntherland is that they are the asymptotic momenta after letting a quantum 
gas expand into a larger volume making it dilute and noninteracting. We exploit this picture to 
make a direct connection to quantities that are accessible in sudden-expansion experiments with 
ultracold quantum gases. By a direct comparison of Bethe-ansatz and time-dependent density matrix 
renormalization group results, we demonstrate that the expansion velocity of a one-dimensional 
Fermi-Hubbard model can be predicted from knowing the distribution of occupied rapidities defined 
by the initial state. Curiously, an approximate Bethe-ansatz solution works well also for the Bose- 
Hubbard model. 


Introduction. Some often emphasized aspects of ex¬ 
periments with ultracold quantum gases as compared to 
condensed matter systems are the high degree of tun- 
ability of dimensionality and interaction strengths and 
the possibility to obtain clean realizations of standard 
many-body Hamiltonians such as the Fermi- and Bose- 
Hubbard model [1-4], albeit usually with inhomogeneous 
density profiles [5]. Moreover, ultracold quantum gases 
are complementary to condensed-matter systems with re¬ 
spect to the typically accessible observables. Primarily, 
these are in-situ density profiles or density profiles after 
time-of-flight experiments (possibly combined with other 
perturbations of the system or quenches of parameters). 
A time-of-flight measurement, achieved by turning off all 
potentials, gives access to either the in-situ momentum 
distribution function or, if the system was initially in 
a lattice and the removal of the lattice occurred suffi¬ 
ciently adiabatically, the quasi-momentum distribution. 
This relies on various assumptions [6, 7], including that 
the atoms do not experience interactions during the time- 
of-flight measurement. 

This latter assumption may either not always be valid 
or one may actually deliberately be interested in eluci¬ 
dating the very effect of interactions, which thus has a 
very different purpose from a time-of-flight experiment. 
Here, we address two questions, namely, first, whether 
the density profiles during the expansion contain any in¬ 
formation about the initial state and, second, what is 
the form of the asymptotic momentum distribution func¬ 
tion. Theoretically, in higher dimensions, the analysis is 
mostly based on time-dependent mean-field schemes [8] 
(see [9, 10] for more recent methodological developments, 
though). In one dimension, for situations that either 
permit a scaling solution (see, e.g., [11, 12]) or are de¬ 
scribed by integrable Hamiltonians such as the Lieb- 
Liniger model of interacting bosons in the continuum, 
these questions can in principle be answered exactly [13- 
19]. We are particularly interested in an integrable lat¬ 
tice model with repulsive onsite interactions, the Fermi- 


Hubbard model (FHM) [20] 

L-l L 
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where L is the number of lattice sites, Cio- is a fermion 
annihilation operator, and Ui^ = cJ^Cicr- In the context of 
bosons, the simplest, completely understood example are 
hard-core bosons (HCBs, the lattice version of the Tonks- 
Girardeau gas [21, 22]), which can be obtained from the 
Bose-Hubbard model (BHM) 

H = -JY^ (al+ifli + h.c.) + 1) (2) 

by sending U/J —t oo or requiring (a|^)^ = 0. Here, 
Oi is a boson annihilation operator and Ui = ajoi. For 
HCBs, the physical quasimomentum distribution func¬ 
tion (quasi-MDF) n/, (defined as the Fourier transform 
of one-particle correlations {a\aj)) after a long expan¬ 
sion time becomes identical to the conserved set 77^ of 
the spinless, noninteracting fermions that HCBs can be 
mapped to [23]. This process goes under the name of 
dynamical fermionization [18, 24] and generalizations ap¬ 
ply to the (integrable) repulsive Lieb-Liniger gas [17, 19] 
and even to the (nonintegrable) BHM with U/J < oo 
[25]. Density profiles {ni{t)) undergo a ballistic expan¬ 
sion for HCBs in a ID lattice, which was observed in 
experiments [26] . The ballistic expansion manifests itself 
in a linear increase R{t) = v^t of the radius defined as 

i 

where N denotes the number of particles and zg = (A -b 
l)/2. The expansion velocity of HCBs is related to the 
conserved 77^ and thus also to the asymptotic form of the 
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physical nk{t = oo) = via [25-27] 

^ ( 4 ) 

k 

where Vk = det/dk are the group velocities of nonin¬ 
teracting particles with a tight-binding dispersion Ck = 
—2Jcos(fc) (the lattice spacing is set to unity). 

In this work, we aim at the generalization of these ob¬ 
servations to Bethe-ansatz (BA) integrable lattice mod¬ 
els with repulsive interactions that do not map onto 
noninteracting particles. Following Sutherland [28], for 
such systems, quasimomenta are replaced by so-called 
rapidities, which have the interpretation that they be¬ 
come physical momenta in the asymptotic regime of an 
expansion. This happens once particles have spatially 
rearranged themselves according to increasing momenta 
and thus stop crossing each other as they continue to ex¬ 
pand. If we then define a distribution of rapidities 
K defined by the initial condition, then our hypothesis 
is that the asymptotic physical momentum distribution 
function nk{t —> oo) becomes equal to the conserved 71^ 
(assuming, for simplicity, real k) 

nfe(too) = . (5) 

As a consequence, since in the asymptotic regime the 
expansion is expected to be ballistic because diluteness 
suppresses any scattering, we expect that the asymptotic 
expansion velocity can be written as 

v'^{t = oo) = (6) 

K 

The main result of our work is that Eq. (6) indeed 
holds for the FHM with repulsive interactions and ini¬ 
tial densities ng = N/Lq smaller or equal than one, ex¬ 
panding from the correlated ground state within a box of 
size Lq (the regime of initial filling larger than one was 
studied numerically in Refs. [29, 30]). Our results are 
based on a comparison of a BA calculation of with 
numerical results for the density profiles obtained from 
time-dependent density matrix renormalization group 
(tDMRG) calculations [31-33]. This implies that a mea¬ 
surement of density profiles, accessible in quantum-gas 
experiments, gives access to the very abstract concept 
of rapidities of an integrable quantum model, which are 
very important for actually carrying out calculations, but 
which are usually hidden. Of course, is just a single 
number and contains only partial information about the 
full rapidity distribution. 

A possible obstacle could be that the times needed to 
reach the asymptotic regime are not accessible to either 
experiments or tDMRG. This is true for the quasi-MDF 
[34, 35], which, using tDMRG, we are only able to ob¬ 
tain for A^ = 2,4 particles in the long-time limit. The 
expansion, however, turns out to be ballistic to a good 
approximation (i.e., R oc t [36]) for the FHM under the 
aforementioned conditions [27] . Such a behavior implies 
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FIG. 1. (Color online) Sudden expansion in the FHM. (al)- 
(a3) Density distribution {uiit)) and (bl)-(b3) renormalized 
quasi-MDF nk{t) at (al), (bl) t = 0, (a2), (b2) tJ = 2, and 
(a3), (b3) tJ = 20 [[/ = 8J, N = 4, Lq = N, we set h = 1], 
(b4) Results for the BHM and the same model parameters. 
Solid lines are tDMRG results, dashed lines in (b3) and (b4) 
show the corresponding Fermi-Dirac function, Eq. (8). In the 
figures, all quantities are expressed in dimensionless units. 


that Vt becomes time independent very rapidly and it 
can thus be extracted already from short-time dynamics, 
long before Uk has converged to its asymptotic regime. 
Thus, experiments do not need to reach the asymptotic 
regime. 

An example, in which Uk nevertheless becomes station¬ 
ary very fast, is the spin-imbalanced FHM with attrac¬ 
tive interactions [35, 37], where a quantum-distillation 
process [29, 38, 39] ensures a fast separation of pairs and 
excess fermions [35]. In that case, the generalization of 
Eq. (5) to both real and complex rapidities (the latter 
present because of the bound states in the ground state 
of the attractive FHM) seems to hold, based on a compar¬ 
ison of tDMRG and BA calculations [35]. Interestingly, 
we will show here that even for the nonintegrable BHM, 
one can exploit a BA approach along the lines of [40] to 
define which via Eq. (6) leads to a good agreement 
with tDMRG data from [25]. 

Asymptotic form of the quasi-MDF at half filling. We 
begin by describing the overall time evolution of densi¬ 
ties and the quasi-MDF, and for the latter, we propose 
a simple expression for its asymptotic form for the ex¬ 
ample of initial states with half filling. Figure 1 shows 
typical results for the FHM at U/J = 8 obtained with 
tDMRG for the density profiles {ni{t)) and the quasi- 
MDF nk(t). We calculated the observables in the initial 
state [Figs. l(al) and l(bl)], in the transient regime of 
the expansion [Figs. I(a2) and l(b2)] and in the asymp¬ 
totic regime, where approaches its stationary form 
[Figs. I(a3) and l(b3)]. The transient regime is charac¬ 
terized by peaks in the quasi-MDF at fc = ±7r/2 [34]. 
Similar transient phenomena for bosons were studied in 
[7, 25, 41, 42]. 
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In the long-time limit, the quasi-MDF of a gas that has 
expanded from a Mott insulator approaches a particle- 
hole symmetric form, both for fermions and bosons [25]. 
For the BHM, this can be viewed as a generalized dy¬ 
namical fermionization [25] , similar to integrable bosonic 
models [18, 19, 24, 43]. The density profile at the longest 
times reached in the simulations is practically flat, except 
for the propagating wavefronts. Therefore, the gas can 
be well approximated by assuming both diluteness and a 
homogeneous density. We find that the final quasi-MDF 
approaches a simple Fermi-Dirac distribution 


where temperature T = 1//3 and chemical potential /r are 
determined to match the energy E, which is conserved 
during the expansion, and the particle number N of the 
strongly correlated system [27]. This effective noninter¬ 
acting gas, containing the same number of particles, can 
be viewed as having originated from the same box. 

For large U/J, the total energy corresponds to rela¬ 
tively high temperatures of the free fermions and, in ad¬ 
dition, /r = 0 at half filling. This simplifies the expression 
of Eq. (7) since the only parameter that determines the 
quasi-MDF of the free particles fk is the energy density 
e = \E\/N. Expressing the quasi-MDE up to 0{e^), we 
get 


fk = 


N 

~L 


^1 -I- ecosfc — 




Note that the first two terms are similar to the results 
of [10]. We use Eq. (8) to compare to the tDMRG re¬ 
sult at fV = 4 for the EHM (e = —0.279J) and the BHM 
(e = —0.364J) in Figs. I(b3) and l(b4), respectively. The 
numerical data away from k = ±7r/2 agree very well with 
the free-fermion reference system. The deviation between 
the two curves can be quantihed by looking at the differ¬ 
ence in average velocities, Auav(i) = 'Cav(i) — , where 

Uav = V2J is the equilibrium average velocity of the 
free-fermion system and Uav(t)^ = (1/fV) 

We find that Avav(t) goes to zero at asymptotic times 
with a power-law dependence. 

While the interpretation of the asymptotic properties 
in terms of a thermal state is very intuitive, it is not 
based on rigorous arguments, unlike the BA approach 
that we will detail next. Nevertheless, the measurement 
of the MDF in an experiment would be very interesting 
and would obtain the entire rapidity distribution. Using 
few particles (as was studied in a recent quantum-walk 
experiment of bosons [44] ) leads to a faster convergence 
towards the stationary form. 

Bethe-ansatz based approach. We turn now to the 
problem of determining the asymptotic expansion veloc¬ 
ity from Eq. (6). Generally speaking, we are attempting 
to predict the asymptotic form of observables from in- 
tegrability in a specific quantum quench problem, which 
is a question that is currently being studied for many 


other quenches, driving methodological advances in the 
framework of the Bethe ansatz (see, e.g., [19, 45, 46]). 
We first need to determine the distribution n^, which is 
conserved for t > 0. The main technical complication 
is that the wavefunction needs to be expanded in the 
postquench eigenstates (after quenching the trapping po¬ 
tential to zero) of the integrable homogeneous FHM in 
an infinitely large lattice. This problem is notoriously 
difficult, and we therefore compute the discrete set of 
prequench rapidities with respect to the initial box (i.e., 
for t < 0, what we denote as ’’rawBA”) and then use 
single-particle projection techniques to approximate 
for t > 0 (’’projBA”). The calculation of the discrete set 
of rapidities for particles in a box is well defined and can 
be done exactly, albeit numerically (see [36, 47] for de¬ 
tails). Next, as the trapping potential is suddenly turned 
off, the initial distribution nK.{t = 0“) is projected into a 
modified one n^it = 0“'') that is consistent with the new 
size and boundary conditions [48] (or lack thereof [49]) 
of the system. Thus, in principle, one needs to compute 
the overlaps between the initial state of the system and 
a complete basis of Bethe states for the system without 
trap. Each of these Bethe states is in one-to-one corre¬ 
spondence with a rapidity distribution and the overlaps 
give the probability amplitudes for combining those into 
the resultant = n„(t = 0+). 

If one can identify the largest overlaps, those will give 
the dominant contributions, while small overlaps will give 
small corrections that can be left out in an approximate 
calculation. On the one hand, for the repulsive models 
we focus on in here, the initial ground states are char¬ 
acterized by having only real-valued charge (and, for the 
EHM, spin) rapidities. On the other hand, both real and 
complex rapidities (strings) are present in the full spec¬ 
trum of these models, but the overlaps in the case of the 
latter are comparatively smaller and can be ignored in 
a first approximation. Complex rapidities are associated 
with different types of bound states and will most often 



FIG. 2. (Color online) Expansion velocity Vr for the Fermi- 
Hubbard model at U/J = 8 with Lq = 20 and the initial 
density no = 0.5, as a function of the reciprocal time, ob¬ 
tained by tDMRG simulations (solid line). In all tDMRG 
simulations, we used the time step At J = 0.01 and the maxi¬ 
mal discarded weight 10“® (different time steps were used to 
check convergence). The black dashed line is a linear fit in 
the time range l/(tJ) < 0.15. 
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FIG. 3. (Color online) Expansion velocity Vr for the FHM. 
Comparison between BA (filled symbols) and tDMRG (open 
symbols). The BA results are obtained with (projBA) or 
without (rawBA) the rapidity projection techniques, (a) Vr 
versus initial density no [Lg = 20; U/J = 2,8]. The dot- 
dashed and the dashed lines are the exact expressions for U = 
0 and U/J = oo, respectively [see Eq. (9)]. (b) Vr versus the 
initial box size Eg [U/J = 8; ng = 0.5]. Lines are fits to the 
rawBA and the tDMRG data in the range l/Lg < 0.1. 


expand more slowly, so by neglecting their contribution 
we overestimate the expansion velocity of the cloud [50] . 
For an initial density approaching the value of one par¬ 
ticle per lattice site, the system gets closer to a Mott- 
insulating state for which double occupancy is relatively 
suppressed (though still non-zero) and thus the projec¬ 
tion onto bound states is also expected to be relatively 
suppressed. The approximation of ignoring the bound- 
state contributions should thus become better the closer 
the trapped system is to no = 1; consistent with our 
numerical comparisons to be shown below. To further 
simplify the calculation, we consider the discrete set of 
initial rapidities treated individually as one-particle dis¬ 
tributions and then combine the resulting post-quench 
distributions to get the final result (this was already seen 
to correctly capture the leading contributions and give 
good results in other situations [35, 43]). Finally, differ¬ 
ent moments of the distribution are combined accord¬ 
ing to Eq. (6) to yield the asymptotic Vr- 

Comparison with tDMRG. The numerical calculation 
using tDMRG proceeds in a very different way. It per¬ 
forms a unitary time evolution of the many-body wave- 
function (approximated via matrix-product states [51]), 
obviously without explicitly connecting to any nontriv¬ 
ial integrals of motion. The expansion velocity of the 
atom cloud converges rapidly to its asymptotic form such 



FIG. 4. (Golor online) Expansion velocity Vr for the BHM. 
Gomparison of data for Vr versus ng between BA (filled sym¬ 
bols) using rapidity projection (projBA), and tDMRG (open 
symbols) [U/J = 8, Eg = 10]. The dashed line is the exact 
expression for U/J = oo, see Eq. (9). The dot-dashed line is 
the result in the dilute limit Vr/J = 2g/E/{JN). 


that only relatively short times need to be reached in 
the simulations [25, 27]. From the tDMRG simulations, 
one can extract the time-dependent expansion velocity 
by calculating Ur(t) = [R{t +dt/2) — R{t —dt/2)]/dt. For 
the systems considered here, Vr{t) does not change con¬ 
siderably in time since the largest difference may be of 
the order of a few percent. A typical time evolution of 
Ur(t) is shown in Fig. 2 for the Fermi-Hubbard model at 
U/J = 8,no= 0.5 and Lq = 20. After a few tunneling 
times oc 1/J, we observe that Vr{t) approaches its asymp¬ 
totic value as ~ l/{tJ). We then obtain the asymptotic 
value of Ur by applying a linear ht Ur(t) = Ur + a/ (tJ) in 
the time interval l/{tJ) 1. 

Figure 3(a) shows the results for the FHM for two dif¬ 
ferent values of the interaction U/ J = 2,8 and Lq = 20. 
The agreement between the tDMRG and BA calculations 
is generally very good. Both approaches consistently 
show that the maximum values of v^/J for strong inter¬ 
actions are found at initial in-trap densities 0.5 < no < 1 
(at n-o = 1 one has Ur/ J = regardless of the interac¬ 
tion strength [27]). For guidance, we also show the exact 
results for Ur for the noninteracting and the C//J —> oo 
limits, computed by taking Lg —)■ oo [27], 


Ur/J 


sin (fcp) cos (fcp) 


1 - 


kp 


(9) 


Here, the Fermi momentum kp is set by the initial 
density ng. For a noninteracting two-component gas, 
kp = ng7r/2, while for the single-component gas that de¬ 
scribes the charge dynamics in the limit of infinite onsite 
repulsion, kp = ng7r. 

In Fig. 3(b), we compare the BA results obtained with 
or without the rapidity projection techniques [36], with 
the tDMRG results. Remarkably, the finite-size scaling 
with respect to the initial box size Lq shows that all data 
sets approach the same value as Lq increases (see also 
[36] for the BHM case). Thus, the approximations used 
in the BA-based approach become increasingly unimpor¬ 
tant. That the asymptotic and the scaling limits of the 












5 


expansion velocities coincide is one of the main nonintu- 
itive findings of our work. 

The case of the BHM is more delicate and already 
within the tDMRG framework, one needs to make a (con¬ 
trolled) approximation by truncating the size of the local 
Hilbert space by introducing a maximum allowed number 
of bosons per site Nc^t- This, as expected, works better 
the stronger the repulsive interaction is (we set fVcut = 5 
for 17/J = 8). From the point of view of BA, the system 
is known to be nonintegrable. Curiously, a system of BA 
equations exists [52] that yield an approximate solution 
which gets also more and more accurate as the interaction 
strength increases, and eventually becomes exact in the 
hard-core limit [40] . Coincidentally, our scheme based on 
the BA rapidities works also the best for strong interac¬ 
tions. We can thus proceed in a similar way as for the 
FHM and compute the expansion velocity Vr as a func¬ 
tion of the initial in-trap density no of the gas. Those 
results are shown in Fig. 4 for a moderately large value 
of 17/J = 8 showing good agreement, in particular, for 
larger initial densities. In addition, the results are quan¬ 
titatively comparable to those for the fermionic case in 
Fig. 3(a) (this is obvious for U/J = oo, see the dashed 
line in Fig. 4, where the density dynamics of fermions 
and bosons is identical [23]). 

The approximate BA equations for the BHM are also 
reliable in the dilute limit, as they tend to the corre¬ 
sponding equations for the Lieb-Liniger model, where in- 
tegrability is restored. In this limit (see the dot-dashed 
line in Fig. 4), it follows from Eq. (6) and the stan¬ 
dard BA expression for the energy of the system that 
Vi-/J = 2yjE/{JN), where E is the energy of the system 
(as calculated in its prequench state and measured with 
respect to the bottom of the tight-binding band); in pre¬ 
cise agreement with the exact result for the Lieb-Liniger 
model [17]. The reason for the recovery of the exact re¬ 
sult is two-fold: (i) there are no bound states (and thus 
no complex rapidities) in the continuum limit as that 
is a lattice effect, and the single-rapidity approximation 
becomes more accurate in the dilute limit. In addition, 
(ii) E is constant after the quench. These considerations 
also apply to the dilute limit of the FHM, so the same 
relation is expected to hold for a Caudin-Yang gas. In 
all cases studied here, because of the repulsive interac¬ 
tions, the asymptotic free Hamiltonian is fermionic [25]; 
even for bosons, in the Lieb-Liniger case, the underlying 
time dependence is captured by knowing that of an an¬ 


tisymmetric free-fermion-like wavefunction characterized 
by the values of the rapidities [15, 53]. 

Discussion. We showed how sudden expansion experi¬ 
ments, already at short times, give access to information 
about integrals of motion that are usually hidden in the 
structure of the wavefunction. Actual experiments of this 
type using cold-atom setups have already been carried 
out for fermions [54] and bosons [7, 21, 26, 39, 55, 56] 
and more accurate ones for both bosons and fermions 
could be within reach exploiting single-site resolution 
techniques [57, 58] (see [44, 59, 60] for work in this di¬ 
rection). The fact that most experiments use a harmonic 
trap does not change the picture qualitatively, it leads 
to a different rapidity distribution in the initial state 
but after removal of the trap, the system is again in- 
tegrable. While we have substantiated the validity of the 
approximations in the BA calculation by a comparison 
to tDMRC for few particles, one can push the BA cal¬ 
culation of asymptotic quantities to much larger initial 
system sizes or particle numbers [36]. 

We speculate that systems close to an integrable point 
are constrained at short times by the integrals of motion 
of that point and reach asymptotic expansion states that 
reflect them since the gas becomes increasingly more di¬ 
lute as it expands [19, 25]. Therefore, perturbations from 
integrability have no time to generate deviations, similar 
to but more robustly so than within the prethermalization 
scenarios realizable on the same type of systems under 
different experimental conditions [61-65]. This conclu¬ 
sion is corroborated by our results for the nonintegrable 
BHM. The identification and study of asymptotic expan¬ 
sion states seems to be a very fertile ground to explore 
the physics of nonequilibrium systems as they constitute 
a very special case of asymptotic states of effectively non¬ 
interacting systems that nevertheless contain much of the 
information about the character and correlations on the 
parent equilibrium states of the actual interacting sys¬ 
tems. 
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Supplemental material 


SI. BETHE-ANSATZ EQUATIONS 

For reference, we give here the Bethe-ansatz equations 
with open boundary conditions that were used to calcu¬ 
late the results for the comparisons of Figs. 3 and 4. For 
the Fermi-Hubbard model, one has two sets of equations, 
one for the charge rapidities, 

TT (sin Kj +Ar, + i^) (sin Kj - ^ 

rhA (sin Kj -I- (sin Kj - 

(SI) 

(where N is the number of particles and M is the number 
of down-spin fermions in the system), and another for the 
spin rapidities, 


TT (An -f Am -|- ic) (A„ — Am -|- ic) 

(An -f Am ic) (An Am ic) 

m^n 

^ A (An+sinAC^ +if) (An-sinAC^ +if) 
(An + sin Kj — if) (An — sin Kj — if) 


The spin and charge rapidities clearly depend on each 
other, yet only the charge rapidities will enter into Eq. (6) 
of the main text. For the Bose-Hubbard model, there are 
only charge rapidities and thus a single set of N equa¬ 
tions. 


N 

n 

m—l 

m^j 


sin Kj — sin Km — ic sin Kj + sin Km — ic 
sin Kj — sin Km + ic sin Kj + sin Km + ic 




(S3) 


In both cases the interaction parameter is c = U/{2J) 
and the total energy is given by the standard formula 
E = —cos Ati. These equations can be brought 
into polynomial form [47] and then solved using stan¬ 
dard root-finding procedures. That process can be done 
more efficiently by solving the equations sequentially and 
following an iterative scheme. Note that in the thermo¬ 
dynamic limit, one could have also assumed the system 
to be in an equilibrium thermal state and determine the 
rapidities by solving the corresponding Thermodynamic 
Bethe-Ansatz equations. 


S2. FINITE-SIZE SCALING 

Given the size limitations of the calculations we do us¬ 
ing tDMRG and, to a different extent, also Bethe ansatz, 



FIG. SI. (Color online) Expansion velocity v,- for the Bose- 
Hubbard model with initial density no = 1 as a function of 
the reciprocal size of the initial system (1/Lo). Comparison 
of calculations based on Bethe ansatz without and with pro¬ 
jection of the rapidities (filled symbols) and using tDMRG 
(open symbols) for a system with U/J = 20. See Fig. 3(b) in 
the main text for an example of the finite-size dependence in 
the FHM case. 


it is interesting to perform a finite-size scaling analysis 
of the results for the expansion velocities. For concrete¬ 
ness, we consider the Bose-Hubbard model and focus on 
the case of an initially half-filled system. As mentioned in 
the main text, for the hq = 1 case, we expect v^/J = \f2 
from numerical studies [25, 27] and the analysis of lim¬ 
iting cases ([/ = 0, U/J = oo) [25, 27]. In Fig. SI we 
show three different calculations, all three converging to 
the quoted expected result. 

The three ways of calculating scale very differently 
with Lq. On the one hand, tDMRG immediately gives an 
accurate result for very small systems that provides up to 
four significant digits of precision; but for larger systems 
the cost of the time evolution needed to approach the bal¬ 
listic expansion regime becomes larger and the error bars 
increase degrading the accuracy of the results. The prac¬ 
tical size limit in this case is about Lq « 20. On the other 
hand, a naive Bethe-ansatz calculation gives results that, 
for small systems, are just order of magnitude estimates; 
but we understand that the reason is that one needs to 
project the rapidities into those of the infinite-size sys¬ 
tem. Performing this projection in the approximate way 
described in the main text yields up to three or four digits 
of precision for small systems (comparable to tDMRG), 
but as the system size grows, that precision degrades to 
only two digits, or a relative error of a few percent. This 
can be ascribed to the growing weight of the bound states 
in the expanding state of the system. Fortunately, within 
the Bethe-ansatz scheme we are using, we have different 
limitations on the size of the systems we can study and 
are able to consider much larger systems (in here we show 
up to Lq = 1000, but even larger numbers of lattice sites 
are possible). The Bethe-ansatz results, with or without 
projection of the rapidities, scale towards the same result. 
In fact, the Lq = 10 case turns out to be the worst-case 
scenario in terms of accuracy and for larger systems we 
again achieve three or four digits of precision. A linear 
regression analysis in terms of 1/Lq yields five or more 
digits of precision (depending on how one restricts the 
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data to larger system sizes). 

The initial state of the Bose-Hubbard model with uq = 
1 is also particularly interesting because we can use the 
properties of the Bethe-ansatz equations to derive the 
exact scaling expression: Vr/J = a/2\/1 -I- 1/Lq. All the 
corresponding data in Fig. SI (green circles) agree with 


this expression to the level of machine precision. From 
there it follows that v^/J « V2+ {l/V^)/Lo -1-C>(1/Lq), 
and we verified that not only the ordinate but also the 
slope of our linear regression are in good agreements. 



